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Dynamic mean field theory is applied to the problem of forest fires. The starting point 
C*") ■ is the Monte Carlo simulation in a lattice of million cells. The statistics of the clusters 

04 ' is obtained by means of the Hoshen-Kopelman algorithm. We get the map p n — > p n +i, 

where p„ is the probability of finding a tree in a cell, and n is the discrete time. We 
demonstrate that the time evolution of p is chaotic. The arguments are provided by 
the calculation of the bifurcation diagram and the Lyapunov exponent. The bifurcation 
diagram reveals several windows of stability, including periodic orbits of length three, five 
and seven. For smaller lattices, the results of the iteration are in qualitative agreement 
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1 Introduction 
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S I It would be of obvious interest if one could predict a forest fire. An attractive aim would be 

i-^j ■ also to evaluate the risk of the forest fires in a given year, having the statistics of the fires in 

the past. Here we argue that the time dependence of the forest area burned in a given year 
is inherently chaotic. This means that such predictions are not possible in longer time scale. 
The method used is the mean field equation based on an appropriate cellular automaton. 

Cellular automata (CA) is a modern tool for qualitative simulations of numerous prob- 
lems in statistical mechanics and theory of dynamical systems [0, ^) . In general, the rules of 
CA are probabilistic — they may depend on a random variable. When this dependence is 
absent, the automaton is termed deterministic. These automata are most exciting but least 
predictable. The behavior of a probabilistic cellular automaton can be predicted successfully 
H by means of a nonlinear mean field equation. The solution is the time-dependent probabil- 
ity distribution of cell states. The approach is termed "mean field" , because the correlations 
between the states of neighboring cells are lost within one time step. However, they are lost 
also in the original probabilistic automaton after several time steps. Similar technique can 
be applied also to investigate deterministic CA ||, but in this case the results of the mean 
field approximation are less reliable: the rules of an original deterministic automaton can 
preserve the correlations forever. 

Our aim here is to apply this mean field technique to the problem of the forest fires. 
The list of references of the problem is vast and reveals its numerous aspects, from purely 
computational to inherently biological ones Q. In particular, some controversy has been 
mentioned Q on the applicability of the percolation theory Q , and on the size dependence 
of results of possible experiments. The advantage of the mean field approximation is that the 
finite-size fluctuations are eliminated. Then, our results can be seen as a limit case both for 
simulations and experiments. Below we demonstrate that the critical concentration for the 
percolation problem is fairly reproduced by the forest fire simulation, if the lattice is large 
enough. 
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The approach can be formulated as follows: The model forest is a two-dimensional regular 
lattice. Each cell at the lattice may be occupied by a tree with probability p. During a year 
(one time step), trees appear at empty cells with a probability rp(l — p), where p is the 
previous concentration of trees per cell, and < r < 1 is a model control parameter. This 
nonlinear character of the probability reflects the obvious tendency that a new tree appears 
less likely if there is almost no free space around the cell. This kind of dependence is known as 
the Verhulst term All trees, those existing previously and those just grown, form clusters 
of different sizes. A cluster is a set of trees which are connected by nearest-neighbor bonds in 
the von Neumann neighborhood. Now, a probability that a cluster is ignited is proportional 
to its size. Then, a fraction of the forest burned each summer is just the mean square of the 
cluster size. We will show that below the percolation threshold p c , this fraction is negligible. 
Above this concentration, most of the forest disappears — only small clusters survive. 

Summarizing our method, the time evolution equation for the concentration of trees in 
the forest is a superposition of two maps, which describe a sequence of growing ("spring") 
and fire ( "summer" ) 

Pn+l = W{p n +rp n (l -Pn)), (la) 

where the map W(p) is 

W{p) = P - A( P ) (lb) 

and n enumerates years, i.e. the time steps. The function A{p) is the burned fraction of the 
total area of the forest, i.e. the mean square of the cluster size. This fraction is evaluated by 
means of the Hoshen-Kopelman algorithm R] . 

The purpose of this paper is to demonstrate, that the map p n — ► p„+i leads to chaotic 
behavior. We note that such a conclusion is not entirely new. In 1990, Chen, Bak and 
Jensen [0 proposed a deterministic coupled map lattice [H for a model description of forest 
fires. In this model, fire was introduced as a continuous process. The Lyapunov exponent 
calculated there described the time evolution of the Hamming distance between two nearby 
trajectories. Initial distance was set by a random perturbation, added to each site at time to. 
As the distance was found to increase with the power law, the value of the Lyapunov exponent 
was zero. Later on, this model was generalized by Socolar, Grinstein and Jayaprakash [jioj] 
by adding a parameter, which allowed to get a positive Lyapunov exponent. 

Both the above mentioned approach and this work have in common that they are com- 
putational rather than oriented to experiment. However, within this computational stream 
there are serious differences. The algorithm of Chen et al. || is deterministic. The forest 
fire is continuous; a tree burns spontaneously as soon as it is large enough. Our map can be 
treated as an extrapolation of the output of a probabilistic automaton to the area of infinite 
size. Then, a sequential iteration is applied, and the season of growing and the season of 
fires are treated as two successive parts of one step of the iteration. From Ref. Q we know 
that there are many uncontrolled quantities relevant to the problem: moisture, humidity, 
landscape, etc. As it is impossible to take them into account within a simple scheme, our 
probabilistic approach seems to be closer to reality. Last but not least, our approach gives a 
sequence of yearly burned areas which, although expressed in arbitrary units, can be directly 
compared to available statistical data. 

The paper is organized as follows: In Section || we describe the way how the Hoshen- 
Kopelman is used here. There, two different kinds of statistics are assigned to casual (light- 
nings) and intentional ("human factor") forest fires. In Section^ the results on the Lyapunov 
exponent and the bifurcation diagram are given. A remark is also added on the fire size dis- 
tribution. The comparison of our results to the statistical data on forest fires in Canada in 
years 1970-2000 Jl2| is included in Section ^. The last section is devoted to discussion. 



2 Numerical approach 

With the Hoshen-Kopelman algorithm |t| we are able to label all occupied sites on a L x L 
large square lattice in such a way that the sites with the same label belong to the same cluster 
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Figure 1 : Are black sites at the line x recognized as a part of the cluster? Algorithm goes 
from left to right, from top to bottom. 



and different labels are associated to different clusters. In this way we get the cluster size 
distribution, which allows us to obtain the function A(p) and, in the next step, to formulate 
the mean field map given in Eq. ([!]). 

The beauty of the Hoshcn-Kopclman algorithm for the cluster characterization is that it 
goes through the lattice only once, and stores only the current line of the length L. Thus 
it is particularly useful for the bus-bar percolation problem investigation and it allows to 
check whether a given occupied site at distance x from an arbitrary chosen (usually the 
first and full occupied) line is still connected to this line through the occupied sites in lines 
2,3,... ,x-2,x-l pi, ||. 



However, to construct the map (lb) it is necessary to evaluate the mean cluster size (or 
its average mass) L 2 A(p). This requires an information on all bonds to sites at the distance 
x from the first line also through the occupied sites at the lines x + 1, x + 2, . . . , L — 1, L, as 
it is presented on Fig. |. This in turn requires storing of the whole lattice for all the time 
of simulation and passing through the lattice twice — what wastes computer memory and 
machine time, but ensures proper labeling of the clusters. 

The simulation is carried out on 10 3 x 10 3 square lattice, with a fraction of p + rp{\ — p) 
occupied sites for different values of the initial trees concentration < po < 1. An average 
concentration of trees in a "green" (but soon "red" and consequently "black" ) cluster 

i 

may be evaluated with two kinds of weights 

w'i = Si/L 2 and w" = Si/ ^ Sj, 

i 

which describe the probability of finding a cluster of size Si on the lattice and the probability 
of choosing an i-th cluster among all other clusters, respectively. These statistics can be 
interpreted as reflecting the probability of a casual fire started from a lightning, which can 
strike at an empty cell (u/), or of a fire started intentionally at a randomly selected tree 
(w"). The obtained map ([!]) for each p is the average over a hundred of independent lattice 
realizations. 

Below we concentrate on the case of the fires of natural origin. The results for the other 
case will be only listed. The details of the latter calculation will be discussed more thoroughly 
elsewhere. 



3 Results 
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Figure 2: Map for model with average burned area A{p) evaluated for different weights w' 
(natural) and w" (provoked). Before percolation threshold burned areas are negligible and 
the map is given by argument of function W in Eq. JTa]). Also map for recipe (|^) is included. 



3.1 The map 

The map p n — > p n +i is given in Fig. 0. Two plots refer to two kinds of statistics, described 
in Section ^. We see that the plots are the same below the maximum of the curve, which 
is near the percolation limit for the square lattice p c = 0.59273 ||, EH). Actually, below this 
limit both plots are identical to the curve p + rp(l — p) from Eq. (|la|). This means, that for 
p < p c the fires are negligible. Above the percolation limit, the curves remain similar for p 
close to p c , but for higher p the curve for the weight w' is remarkably higher, than the curve 
for the weight w" . It is clear that intentional fires (w") are more devastating, because in this 
case it is sure that a tree is ignited. 



3.2 The bifurcation diagram 

The rest of the results refer to the case of casual fires, i.e. the weight w' . In Fig. ||we show 
the bifurcation diagram. This diagram is obtained with the map p n — > p n +i, supplemented 
by straight lines joining neighboring 10 3 numerical points. We can see the windows of stable 
cycles of lengths four, eig ht, fiv e, ten, six, twelve, seven and eight again (Fig. 3(a)| ), five, 
three, five and two in Fig. 3(b) , when the parameter r decreases. The presence of the cycle 



three proves that, by means of the Sharkovskii theorem, all possible lengths of the limit cycles 
are present 16 1. 

A question arises, if our bifurcation diagram is not an artificial consequence of the finite 
grid of the map and/or the finiteness of the lattice. To check this point, we calculated also 
the bifurcation diagram for a piecewisely analytical map G(H(p)), where 



H(p) = p + rp(l -p), 



G(p) = 



0.1 exp[100(0.6 - p)} + OA^/T^p 



p > 0.6, 
p < 0.6. 



(2) 
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Figure 3: Continued from page|^. 



This map is also shown in Fig. |[ It is designed not as to fit the plot obtained within the 
Hoshcn-Kopclman algorithm, but rather to investigate the co nseq uences of the particular 
shape of this plot. The obtained bifurcation diagram (Fig. |3(c)| ) shows the windows of 
stability, with the width which increases with the parameter r. The lengths of the cycles 
decreases with r in the same way as presented in Fig. 3(a) . We deduce that this part of the 
bifurcation diagram is generic. 

To each cycle, an admissible word can be assigned according to the rules of symbolic 
dynamics. We have only to check the values of p in the sequence of a given limit cycle. The 
rule is as follows: The value of p preceding the largest p (close to p c ) in the cycle is labeled 
C. Let us call this value p* . The next value (close to p c ) is labeled R. Any other p is labeled 
R if it is larger than p* , and it is labeled L if its value is smaller than p* . We use this 
rule to assign symbolic words to each cycle in Fig. 3(a), and we compared their sequence 



to the standard one for the unimodal maps [fL7| . The sequences of the cycles presented in 
Figs. [3(a)] and gb| are respectively RL 2 , RLrRL 2 , RL 3 , RL^RL 3 , RL A , RL 5 RL 4 , RL 5 , 
RL 6 and RL 2 R, RL, RLR 2 , R with decreasing parameter r. Note that C is omitted for 
convenience at the beginning of each word. We have found that these sequences indicate 
opposite directions of the increase of the parameter r. According to the results discussed in 
the preceding paragraph, the left part of the bifurcation diagram (Fig. |3(b)[ ) is supposed to 
vanish in the thermodynamic limit. Still, the conclusion on all possible lengths of the limit 
cycles holds, as it follows from the presence of cycle five (all lengths but three) and cycle seven 
(all but three and five) p| . For the bifurcation diagram in Fig. 3(c) the above prescription 
gives words RL n . However, in this case the maximum of the map does not appear in the 
limit cycle, except at the bifurcation points where the cycle length changes. In this sense the 
letter C should not be used. 

We have checked that the bifurcation diagram for weigths w" seems to be a homogeneous 
spot. Still, we have found the cycle of length three by the symbolic dynamics method fTlf| . 
The details will be given elsewhere. 



6 



0.4 



0.3 



0.2 



0.1 



<< - 



-0.1 



-0.2 



-0.3 



-0.4 



w' (natural), p =0.5, e=10 15 ...10" 10 , N 



pts" 



=100 




t 



J? 



+++ 



++* + v+v+ 



+ + 

+ + * 
+♦+ i + 



+ + + +++ ++ + + +^lf% + + * 
+ + + +* \% + + * 



++ 



0.2 



0.4 



0.6 



0.8 



Figure 4: Positive Lyapunov exponent for different model control parameter r. Results of 
evaluating A are averaged over few differences e between two values of initial concentration 
Po- The linear fit to the first hundred points of the Fig. 5(b) were used. 



3.3 The Lyapunov exponent 

In Fig. |] we show the largest Lyapunov exponent as dependent on the parameter r. Only 
positive values are shown, and the spaces between data agree with the windows of stability in 
the bifurcation diagram. We have found that the negative values of the Lyapunov exponent 
are much more difficult to be calculated when we have a limit cycle and not a fixed point. 
In this case, the time dependence of the difference between trajectories in not monotonic. A 
trajectory can be trapped by a limit cycle with a shifted phase, which leads to a fixing of the 
difference between trajectories even in the case of a stable limit cycle. Moreover, sometimes 
this trapping can be only temporal. As the result, in the presence of a stable limit cycle 
the difference between trajectories can depend on the initial point. These complexities are 
revealed in Fig. ||. 

The Lyapunov exponent for weights w" seems to be positive in almost the whole range 
of the parameter r. This confirms that the windows of stability are much more narrow. 



3.4 Fire size distribution 

We have calculated the histogram of the size of fires, given by the function A(p n ) = A n . 
Typical results for the map calculated by Eq. ([!]) for the weights w' are given in Fig. |6(b)| . 
As we see, the "fires of size zero" are most frequent. However, larger fires can be more likely 
that smaller ones. This kind of curves do not show any self-similarity. 

A serious difference between the results for w' and w" is found for the histogram of the 
forest fires. An example for r = 0.3 is presented in Fig. 0. For large fires, the curve seems to 
follow the exponential dependence. To capture this tendency, we have fitted the map for w" 
by a curve similar to Eq. (^), where the function G(p) is chosen as 



G{p) = Pc exp(-12.1735Vp :r ^) (3) 
if p > p c — 0.59237. The obtained histogram is added to Fig. [?]. Note that the latter curve 
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Figure 5: Times dependence of the difference between two trajectories which start from po 
and po + e. Here e = 10~ 15 and (a) r = 0.5, (b) r = 1.0. 
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Figure 6: Histogram of burned area (a) in Canada, 1970-2000 [^2| rescaled to the largest fire 
in a given territory, and (b) obtained by means of computer simulation. 
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Figure 7: Histogram of burned area for weights w" . 



is monotonic except the point at A = 0, i.e. the events without a fire. As we see, small finite 
fires are neglected by the above analytic approximation. 



4 Comparison with experimental data 

Our map p n — > p n +i is obtained for the lattice of million sites, as an average over hundred 
realizations. The accessible statistical data p2) relate to areas larger than million trees each, 
but to thirty "iterations" for one initial state of the forest. On the other hand, several 
uncontrolled factors (moisture, etc.) cannot be taken into account within our simplified 
approach. We expect that these agents lead to some noise, which is absent in our deterministic 
map. This noise can be reproduced here by decreasing the model lattice. That is why we are 
going to compare the experimental data |l2| with the results of our simulation performed for 
a small lattice. 

We have checked that the time dependence of the size of the fires of both experiment and 
calculations show the same kind of oscillations. Examples are given in Fig. [|. In Fig. 8(a) 



the data are presented on the forest fire in Yukon, Canada, in years 1971-1998 |12|. In Fig. 
gbjthe 

same data are shown for the simulation performed on the lattice of 50 x 50 cells. As 
we see, the qualitative features of the data are reflected by the simulation. 

We also checked that the data on other provinces in Canada in years 1970-2000 |j3 do 
not differ in any essential way. In Fig. |(a)] we show a collective histogram for experimental 



data. For each province or territory, the fires are rescaled to the maximal fire in this region 
within the years 1970-2000. The statistics does not allow to speculate on the size distribution 
of large fires. However, it is obvious that very small fires are most frequent. This particular 
effect is reproduced in our numerical results (Fig. |6(b)| and . 

5 Discussion 



The bifurcation diagram for weights w' in Figs. 3(a) - 3(b) reveals an unusual structure, which 



is a conglomerate of two diagrams. On the right side of the plot we observe some period 
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Figure 8: Forest fire statistic for (a) real and (b) computer experiment. 
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doubling when the parameter r decreases. On the left, clear period doublings are observed 
from a fixed point to a limit cycle of length two, and then to a cycle of length four, when the 
parameter r increases. This left part of the bifurcation diagram is due to the size effect of 
finite lattice. On the other hand, the map for weights w" constructed in the same way does 
not reveal wide windows of stability. This means that the period doublings and the windows 
are sensitive to the shape of the map on the right side of the percolation threshold. We could 
only prove that the cycles of all lengths arc present there, by means of the presence of the 
cycle of period three. The map is similar to the tent map, where no windows are observed. 

Concluding this point, the bifurcation diagram reveals that the map is not unimodal, what 
is not a surprise. Its structure may be caused by an intermittent character of the trajectories 
for small r. In this case, if the initial value of p is small, forest fires do not occur in several 
years and the evolution is governed by its parabolic part p + rp(l —p), which remains almost 



constant through a long time. The diagrams in Figs. 3(a) and 3(c) arc similar to the diagram 
obtained in |l8| . It was produced using a discontinuous map, where the time evolution was 
also intermittent. On the other hand, for larger r the flat shape of the map above p c stabilizes 
the limit cycles. We note that the intermittent character of the problem is to some extent 
reflected also by the experimental data, where the majority of fires is of small sizes. 

The histogram in Fig. |^ for weights w" reveals an exponential curve of the number of fires 
of size A in the range of intermediate and large A. The number of fires of size A increases 
with A. This result is in contradiction to what is known from literature on the self-organized 
criticality |l9[| . We interpret it as a particular consequence of our assumption that fires are 
initiated more likely in large clusters. In the approach presented in fl9|| , a tree is ignited 
spontaneously if it is old enough. 

Coming back to the problem of predictability, we note that when a deterministic mean 
field theory in the thermodynamic limit indicates the chaotic behavior of the system, it is 
hard to imagine that any other approach can produce a controllable solution. The exception 
is the the presence of the windows of stability, which allows us to predict the size of a forest 
fire for given values of r. Still, each limit cycle contains large fires, indicated by the letter 
R in the respective word. Obviously, this kind of trajectory is not recommended by forest 
rangers. Concluding, both the bifurcation diagram and the Lyapunov exponent reveal that 
chaos is indeed present in the time evolution of the forest fire. 

It seems also that the mean field approach can be a reference method in investigating other 
problems where probabilistic cellular automata are useful. A list of such problems includes 
immunology |2(i| ] and magnetic systems pl| . To apply a map or a differential equation? The 
answer must be specific for a given problem. 
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